Loading the libraries and installation of the package.

#install.packages("~/Desktop/tutorial/space_0.1-1.1.tar.gz", repos = NULL, type = "source")
#install_github(repo = "topherconley/spacemap")
suppressPackageStartupMessages(library(igraph))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(ggpubr))
suppressPackageStartupMessages(library(space))
suppressPackageStartupMessages(library(glasso))
#suppressPackageStartupMessages(library(spacemap))

Run glasso on example data

Estimates a sparse inverse covariance matrix using a lasso (L1) penalty

glasso(s, rho, nobs=NULL, zero=NULL, thr=1.0e-4, maxit=1e4,  approx=FALSE, 
penalize.diagonal=TRUE, start=c("cold","warm"), w.init=NULL,wi.init=NULL, trace=FALSE)

s   : Covariance matrix:p by p matrix (symmetric)

rho : (Non-negative) regularization parameter for lasso. rho=0 means no regularization.

data(spaceSimu)

Y.m<- spaceSimu$Y.data
p<- dim(Y.m)[2] # no. of nodes
n=nrow(Y.m)
s<- var(Y.m)
a<-glasso(s, rho=0.25) ####### change the rho parameter to control the number of edges.
a$w[1:10,1:10]

##                [,1]          [,2]          [,3]          [,4]          [,5]
##  [1,]  1.250000e+00 -7.133126e-03  6.382393e-03  9.811987e-11  2.015679e-04
##  [2,] -7.133126e-03  1.250000e+00 -4.257136e-05 -1.726507e-08 -1.344484e-06
##  [3,]  6.382393e-03 -4.257136e-05  1.250000e+00 -6.239176e-13  3.947731e-02
##  [4,]  9.811987e-11 -1.726507e-08 -6.239176e-13  1.250000e+00 -4.957146e-14
##  [5,]  2.015679e-04 -1.344484e-06  3.947731e-02 -4.957146e-14  1.250000e+00
##  [6,]  9.372103e-03 -5.348190e-05  4.785325e-05  6.607423e-13  1.511298e-06
##  [7,]  1.042896e-03 -6.956248e-06  2.042524e-01 -2.564615e-13  2.415964e-01
##  [8,]  1.995857e-05 -1.138977e-07  1.019207e-07  3.229361e-14  3.219189e-09
##  [9,] -6.325406e-04  1.108456e-01 -3.775078e-06 -1.531005e-09 -1.192241e-07
## [10,]  3.688676e-03 -2.191352e-05  1.926529e-05  2.938483e-13  6.084340e-07
##                [,6]          [,7]          [,8]          [,9]         [,10]
##  [1,]  9.372103e-03  1.042896e-03  1.995857e-05 -6.325406e-04  3.688676e-03
##  [2,] -5.348190e-05 -6.956248e-06 -1.138977e-07  1.108456e-01 -2.191352e-05
##  [3,]  4.785325e-05  2.042524e-01  1.019207e-07 -3.775078e-06  1.926529e-05
##  [4,]  6.607423e-13 -2.564615e-13  3.229361e-14 -1.531005e-09  2.938483e-13
##  [5,]  1.511298e-06  2.415964e-01  3.219189e-09 -1.192241e-07  6.084340e-07
##  [6,]  1.250000e+00  7.819330e-06  1.496433e-07 -4.742587e-06  2.765652e-05
##  [7,]  7.819330e-06  1.250000e+00  1.665582e-08 -6.168556e-07  3.147987e-06
##  [8,]  1.496433e-07  1.665582e-08  1.250000e+00 -1.010005e-08  5.889660e-08
##  [9,] -4.742587e-06 -6.168556e-07 -1.010005e-08  1.250000e+00 -1.943214e-06
## [10,]  2.765652e-05  3.147987e-06  5.889660e-08 -1.943214e-06  1.250000e+00

a$wi[1:10,1:10]

##               [,1]         [,2]       [,3]     [,4]       [,5]         [,6]
##  [1,]  1.007181838  0.004684408  0.0000000 0.000000  0.0000000 -0.005998483
##  [2,]  0.004684408  0.849759746  0.0000000 0.000000  0.0000000  0.000000000
##  [3,]  0.000000000  0.000000000  0.9017998 0.000000  0.0000000  0.000000000
##  [4,]  0.000000000  0.000000000  0.0000000 0.800471  0.0000000  0.000000000
##  [5,]  0.000000000  0.000000000  0.0000000 0.000000  0.8310445  0.000000000
##  [6,] -0.005998483  0.000000000  0.0000000 0.000000  0.0000000  0.861787423
##  [7,]  0.000000000  0.000000000 -0.1286478 0.000000 -0.1606219  0.000000000
##  [8,]  0.000000000  0.000000000  0.0000000 0.000000  0.0000000  0.000000000
##  [9,]  0.000000000 -0.071502661  0.0000000 0.000000  0.0000000  0.000000000
## [10,]  0.000000000  0.000000000  0.0000000 0.000000  0.0000000  0.000000000
##             [,7]      [,8]        [,9]     [,10]
##  [1,]  0.0000000 0.0000000  0.00000000 0.0000000
##  [2,]  0.0000000 0.0000000 -0.07150266 0.0000000
##  [3,] -0.1286478 0.0000000  0.00000000 0.0000000
##  [4,]  0.0000000 0.0000000  0.00000000 0.0000000
##  [5,] -0.1606219 0.0000000  0.00000000 0.0000000
##  [6,]  0.0000000 0.0000000  0.00000000 0.0000000
##  [7,]  0.8572223 0.0000000  0.00000000 0.0000000
##  [8,]  0.0000000 0.8304466  0.00000000 0.0000000
##  [9,]  0.0000000 0.0000000  0.80634068 0.0000000
## [10,]  0.0000000 0.0000000  0.00000000 0.8186895

Plot the network with detected edges

fit.adj=abs(a$wi)>1e-6
sum(fit.adj==1)/2                      

## [1] 781

temp=graph.adjacency(adjmatrix=fit.adj, mode="undirected")
#plot(temp, vertex.size=3, vertex.frame.color="white",layout=layout.fruchterman.reingold, vertex.label=NA, edge.color=grey(0.5))

Calculate the degree of the nodes and get the hub nodes

Degree of every node = number of edges this node has with other nodes.
Hub nodes: Nodes connected with large number of other nodes.

temp.degree=apply(fit.adj, 2, sum)

tab<- cbind.data.frame("Nodes"=colnames(Y.m), "degree"=temp.degree)
s<- sort(tab$degree, ind=T, decreasing = T)
tab.sort<- tab[s$ix,]
tab.sort$Nodes<- factor(tab.sort$Nodes, levels=unique(tab.sort$Nodes))
dat<- tab.sort[1:30,]

#ggplot(dat, aes(y=degree, x=Nodes)) +geom_bar(position="stack", stat="identity") + theme_bw() + theme(axis.text.x = element_text(colour="black",size=15,angle=45,hjust=1,vjust=1,face="plain"),axis.title.x = element_text(colour="black",size=25,angle=0,hjust=.5,vjust=.5,face="plain"), axis.text.y = element_text(colour="black",size=15,angle=0,hjust=1,vjust=0,face="plain"),  axis.title.y = element_text(colour="black",size=25,angle=90,hjust=.5,vjust=.5,face="plain"),legend.title = element_text(colour="black",size=15), legend.text=element_text(colour="black",size=15), legend.direction = "horizontal", legend.position = "top") + theme(plot.margin = unit(c(1,1,1,1), "cm"))

Run space on example data using space.neighbor

space.neighbor: A function to estimate partial correlations using the neighborhood selection approach

space.neighbor(Y.m, lam1, lam2=0, iter)

lam1: numeric value. This is the l_1  norm penalty parameeter. 
lam2: numeric value. If not specified, lasso regression is used in the neighborhood selection.

data(spaceSimu)

Y.m<- spaceSimu$Y.data
p<- dim(Y.m)[2] # no. of nodes
n=nrow(Y.m)

alpha=1
l1=1/sqrt(n)*qnorm(1-alpha/(2*p^2))
iter=3

result1=space.neighbor(Y.m, lam1=l1*0.7, lam2=0)
result1$ParCor[1:10,1:10]

##       [,1]      [,2]    [,3] [,4]      [,5] [,6]      [,7] [,8]      [,9] [,10]
##  [1,]    1 0.0000000 0.00000    0 0.0000000    0 0.0000000    0 0.0000000     0
##  [2,]    0 1.0000000 0.00000    0 0.0000000    0 0.0000000    0 0.1049529     0
##  [3,]    0 0.0000000 1.00000    0 0.0000000    0 0.1543100    0 0.0000000     0
##  [4,]    0 0.0000000 0.00000    1 0.0000000    0 0.0000000    0 0.0000000     0
##  [5,]    0 0.0000000 0.00000    0 1.0000000    0 0.2620895    0 0.0000000     0
##  [6,]    0 0.0000000 0.00000    0 0.0000000    1 0.0000000    0 0.0000000     0
##  [7,]    0 0.0000000 0.15431    0 0.2620895    0 1.0000000    0 0.0000000     0
##  [8,]    0 0.0000000 0.00000    0 0.0000000    0 0.0000000    1 0.0000000     0
##  [9,]    0 0.1049529 0.00000    0 0.0000000    0 0.0000000    0 1.0000000     0
## [10,]    0 0.0000000 0.00000    0 0.0000000    0 0.0000000    0 0.0000000     1

Plot the network with detected edges

fit.adj=abs(result1$ParCor)>1e-6
sum(fit.adj==1)/2                      

## [1] 729

rownames(fit.adj) <- colnames(fit.adj) <- colnames(Y.m)
temp=graph.adjacency(adjmatrix=fit.adj, mode="undirected")
# plot(temp, vertex.size=3, vertex.frame.color="white",layout=layout.fruchterman.reingold, vertex.label=NA, edge.color=grey(0.5))

Calculate the degree of the nodes and get the hub nodes

Degree of every node = number of the parent nodes + number of child nodes

temp.degree=apply(fit.adj, 2, sum)

tab<- cbind.data.frame("Nodes"=colnames(Y.m), "degree"=temp.degree)
s<- sort(tab$degree, ind=T, decreasing = T)
tab.sort<- tab[s$ix,]
tab.sort$Nodes<- factor(tab.sort$Nodes, levels=unique(tab.sort$Nodes))
dat<- tab.sort[1:30,]
#ggplot(dat, aes(y=degree, x=Nodes)) +
  # geom_bar(position="stack", stat="identity") + theme_bw() +
  # theme(axis.text.x = element_text(colour="black",size=15,angle=45,hjust=1,vjust=1,face="plain"),axis.title.x = element_text(colour="black",size=25,angle=0,hjust=.5,vjust=.5,face="plain"), axis.text.y = element_text(colour="black",size=15,angle=0,hjust=1,vjust=0,face="plain"),  axis.title.y = element_text(colour="black",size=25,angle=90,hjust=.5,vjust=.5,face="plain"),legend.title = element_text(colour="black",size=15), legend.text=element_text(colour="black",size=15), legend.direction = "horizontal", legend.position = "top") + theme(plot.margin = unit(c(1,1,1,1), "cm"))

Extracting modules from the network

Modules are closely connected markers with similar functional roles/activities

obj<- graph_from_adjacency_matrix(fit.adj, mode = "undirected", weighted = TRUE)

formodule<-cluster_edge_betweenness(obj, weights = NULL, directed = FALSE, edge.betweenness = TRUE, merges = TRUE, bridges = TRUE, modularity = TRUE, membership = TRUE)

module<-communities(formodule)
size<- unlist(lapply(module,length))
mod.s5<- module[size > 5]
mod.s5[1:10]

## $`1`
##  [1] "g 1"   "g 10"  "g 11"  "g 12"  "g 13"  "g 37"  "g 46"  "g 54"  "g 66" 
## [10] "g 99"  "g 441" "g 459" "g 474" "g 480" "g 499"
## 
## $`2`
##  [1] "g 2"  "g 19" "g 30" "g 31" "g 58" "g 59" "g 61" "g 79" "g 82" "g 83"
## [11] "g 84" "g 92"
## 
## $`3`
##  [1] "g 3"   "g 5"   "g 7"   "g 14"  "g 21"  "g 36"  "g 42"  "g 49"  "g 50" 
## [10] "g 52"  "g 63"  "g 75"  "g 77"  "g 85"  "g 90"  "g 94"  "g 353"
## 
## $`4`
##  [1] "g 4"   "g 6"   "g 18"  "g 29"  "g 38"  "g 67"  "g 72"  "g 76"  "g 95" 
## [10] "g 98"  "g 401" "g 408" "g 421" "g 427" "g 445"
## 
## $`5`
##  [1] "g 8"   "g 32"  "g 41"  "g 43"  "g 47"  "g 87"  "g 88"  "g 411" "g 418"
## [10] "g 456" "g 468"
## 
## $`6`
##  [1] "g 9"   "g 48"  "g 93"  "g 305" "g 310" "g 320" "g 332" "g 350" "g 358"
## [10] "g 361" "g 363" "g 364" "g 379" "g 389"
## 
## $`7`
##  [1] "g 15"  "g 20"  "g 25"  "g 35"  "g 39"  "g 53"  "g 60"  "g 62"  "g 89" 
## [10] "g 409" "g 417" "g 424" "g 464" "g 484" "g 491"
## 
## $`9`
##  [1] "g 17"  "g 45"  "g 81"  "g 110" "g 120" "g 143" "g 146" "g 174" "g 175"
## [10] "g 178" "g 185" "g 191" "g 198" "g 431" "g 438" "g 458" "g 471" "g 472"
## [19] "g 475" "g 482" "g 498"
## 
## $`11`
##  [1] "g 23"  "g 28"  "g 109" "g 121" "g 124" "g 135" "g 139" "g 144" "g 145"
## [10] "g 152" "g 155" "g 156" "g 164" "g 167" "g 186" "g 189" "g 194"
## 
## $`15`
##  [1] "g 44"  "g 235" "g 307" "g 311" "g 317" "g 328" "g 335" "g 337" "g 351"
## [10] "g 354" "g 355" "g 367" "g 378" "g 384" "g 391" "g 393"

diag(fit.adj)<- F
mod<- fit.adj[which(rownames(fit.adj) %in% module[[1]]), which(colnames(fit.adj) %in% module[[1]])]
sum(mod)

## [1] 30

dim(mod)

## [1] 15 15

p.m=dim(mod)[1]
obj<- graph_from_adjacency_matrix(mod, mode = "undirected", weighted = TRUE)
# pdf("mod.spacenb.pdf",width = 40, height = 20, useDingbats=T)
# plot.igraph(obj,vertex.size=10, edge.length=5, edge.color="black", vertex.label.cex=2.5,vertex.color=V(obj)$color,layout=layout.fruchterman.reingold)
# dev.off()

Run space on example data using space.joint

space.joint: A function to estimate partial correlations minimizing the joint loss function.

space.joint(Y.m, lam1, lam2=0, iter)

data(spaceSimu)

Y.m<- spaceSimu$Y.data
p<- dim(Y.m)[2] # no. of nodes
n=nrow(Y.m)

alpha=1
l1=1/sqrt(n)*qnorm(1-alpha/(2*p^2))
iter=3

result2=space.joint(spaceSimu$Y.data, lam1=l1*n*1.56, lam2=0, iter=iter)

## [1] "iter=1"
## [1] "iter=2"
## [1] "iter=3"

result2$ParCor[1:10,1:10]

##       [,1]       [,2]      [,3] [,4]      [,5] [,6]      [,7] [,8]       [,9]
##  [1,]    1 0.00000000 0.0000000    0 0.0000000    0 0.0000000    0 0.00000000
##  [2,]    0 1.00000000 0.0000000    0 0.0000000    0 0.0000000    0 0.09927696
##  [3,]    0 0.00000000 1.0000000    0 0.0000000    0 0.1270736    0 0.00000000
##  [4,]    0 0.00000000 0.0000000    1 0.0000000    0 0.0000000    0 0.00000000
##  [5,]    0 0.00000000 0.0000000    0 1.0000000    0 0.2424002    0 0.00000000
##  [6,]    0 0.00000000 0.0000000    0 0.0000000    1 0.0000000    0 0.00000000
##  [7,]    0 0.00000000 0.1270736    0 0.2424002    0 1.0000000    0 0.00000000
##  [8,]    0 0.00000000 0.0000000    0 0.0000000    0 0.0000000    1 0.00000000
##  [9,]    0 0.09927696 0.0000000    0 0.0000000    0 0.0000000    0 1.00000000
## [10,]    0 0.00000000 0.0000000    0 0.0000000    0 0.0000000    0 0.00000000
##       [,10]
##  [1,]     0
##  [2,]     0
##  [3,]     0
##  [4,]     0
##  [5,]     0
##  [6,]     0
##  [7,]     0
##  [8,]     0
##  [9,]     0
## [10,]     1

Plot the network with detected edges

fit.adj=abs(result2$ParCor)>1e-6
sum(fit.adj==1)/2                      

## [1] 731

rownames(fit.adj) <- colnames(fit.adj) <- colnames(Y.m)

temp=graph.adjacency(adjmatrix=fit.adj, mode="undirected")
# plot(temp, vertex.size=3,vertex.frame.color="white",layout=layout.fruchterman.reingold, vertex.label=NA, edge.color=grey(0.5))

Calculate the degree of the nodes and get the hub nodes

Degree of every node = number of the parent nodes + number of child nodes

temp.degree=apply(fit.adj, 2, sum)

tab<- cbind.data.frame("Nodes"=colnames(Y.m), "degree"=temp.degree)
s<- sort(tab$degree, ind=T, decreasing = T)
tab.sort<- tab[s$ix,]
tab.sort$Nodes<- factor(tab.sort$Nodes, levels=unique(tab.sort$Nodes))
dat<- tab.sort[1:30,]
# ggplot(dat, aes(y=degree, x=Nodes)) +
#   geom_bar(position="stack", stat="identity") + theme_bw() +
#   theme(axis.text.x = element_text(colour="black",size=15,angle=45,hjust=1,vjust=1,face="plain"),axis.title.x = element_text(colour="black",size=25,angle=0,hjust=.5,vjust=.5,face="plain"), axis.text.y = element_text(colour="black",size=15,angle=0,hjust=1,vjust=0,face="plain"),  axis.title.y = element_text(colour="black",size=25,angle=90,hjust=.5,vjust=.5,face="plain"),legend.title = element_text(colour="black",size=15), legend.text=element_text(colour="black",size=15), legend.direction = "horizontal", legend.position = "top") + theme(plot.margin = unit(c(1,1,1,1), "cm"))

Extracting modules from the network

Modules are closely connected markers with similar functional roles/activities

obj<- graph_from_adjacency_matrix(fit.adj, mode = "undirected", weighted = TRUE)

formodule<-cluster_edge_betweenness(obj, weights = NULL, directed = FALSE, edge.betweenness = TRUE, merges = TRUE, bridges = TRUE, modularity = TRUE, membership = TRUE)

module<-communities(formodule)
size<- unlist(lapply(module,length))
mod.s5<- module[size > 5]
mod.s5[1:10]

## $`1`
##  [1] "g 1"  "g 8"  "g 10" "g 11" "g 12" "g 13" "g 16" "g 32" "g 33" "g 34"
## [11] "g 37" "g 41" "g 43" "g 46" "g 47" "g 54" "g 66" "g 80" "g 87" "g 88"
## [21] "g 99"
## 
## $`2`
##  [1] "g 2"   "g 9"   "g 17"  "g 19"  "g 30"  "g 31"  "g 45"  "g 48"  "g 56" 
## [10] "g 58"  "g 59"  "g 61"  "g 72"  "g 79"  "g 81"  "g 82"  "g 83"  "g 84" 
## [19] "g 92"  "g 93"  "g 100"
## 
## $`3`
##  [1] "g 3"  "g 5"  "g 7"  "g 21" "g 36" "g 49" "g 50" "g 52" "g 63" "g 90"
## [11] "g 94"
## 
## $`4`
##  [1] "g 4"   "g 29"  "g 67"  "g 75"  "g 76"  "g 85"  "g 95"  "g 98"  "g 349"
## [10] "g 353"
## 
## $`5`
## [1] "g 6"   "g 18"  "g 38"  "g 401" "g 408" "g 421" "g 427" "g 445" "g 478"
## 
## $`7`
## [1] "g 15" "g 20" "g 25" "g 39" "g 53" "g 62"
## 
## $`9`
## [1] "g 23"  "g 28"  "g 121" "g 124" "g 139" "g 144" "g 156" "g 164" "g 189"
## 
## $`14`
##  [1] "g 44"  "g 60"  "g 307" "g 311" "g 324" "g 328" "g 355" "g 378" "g 384"
## [10] "g 391" "g 393"
## 
## $`18`
##  [1] "g 101" "g 105" "g 106" "g 109" "g 116" "g 122" "g 125" "g 128" "g 130"
## [10] "g 132" "g 135" "g 148" "g 152" "g 155" "g 158" "g 160" "g 163" "g 176"
## [19] "g 180" "g 184" "g 195" "g 200"
## 
## $`19`
##  [1] "g 102" "g 108" "g 113" "g 119" "g 127" "g 131" "g 134" "g 149" "g 150"
## [10] "g 154" "g 157" "g 168" "g 182" "g 192" "g 196"

Run space on example data using spacemap

# load("simdata.Rdata")
# p<- dim(sim1$Y)[2] # no. of nodes
# colnames(sim1$Y)<- paste0("Y_",1:p)

# net <- spacemap(Y = sim1$Y, X = sim1$X, lam1 = 70, lam2 = 28.8, lam3 = 12.38)
# net$ParCor[1:5,1:5]
# net$sig.fit[1:10]
# 
# adjnet <- adjacency(net)
# dim(adjnet$yy)
# dim(adjnet$xy)
# 
# adjnet$yy[1:5,1:5]
# adjnet$xy[1:5,1:5]
# adj.yy<- adjnet$yy
# adj.xy<- adjnet$xy
# 
# colnames(adj.yy)<- colnames(sim1$Y)
# 
# obj<- graph_from_adjacency_matrix(adj.yy, mode = "undirected", weighted = TRUE)

# pdf("netyy.pdf",width = 30, height = 30, useDingbats=T)
# plot.igraph(obj,vertex.size=10, edge.length=5, edge.color="black", vertex.label.cex=4,vertex.color=V(obj)$color,layout=layout.fruchterman.reingold)
# dev.off()

Ranking of hub nodes

# adj<- adj2igraph(adj.yy, adj.xy, yinfo = NULL, xinfo = NULL, weighted = NULL)
# 
# r.hub<- rankHub(adj, bdeg = NULL, level = c("x", "y"))
# r.hub[[1]]
# r.hub[[10]]
# r.hub[[15]]

Applying glasso and space on CPTAC-2016 protein data

Running glasso

load("data_retro.RData")  
data<- t(data)
p<- dim(data)[2] # no. of nodes
n=nrow(data)
s<- var(data)
a<-glasso(s, rho=0.05)
a$w[1:10,1:10]

##             [,1]       [,2]      [,3]       [,4]       [,5]         [,6]
##  [1,] 0.08227358 0.00000000 0.0000000 0.00000000 0.00000000 0.000000e+00
##  [2,] 0.00000000 0.09158301 0.0000000 0.00000000 0.00000000 0.000000e+00
##  [3,] 0.00000000 0.00000000 0.1060627 0.00000000 0.00000000 0.000000e+00
##  [4,] 0.00000000 0.00000000 0.0000000 0.08583836 0.00000000 0.000000e+00
##  [5,] 0.00000000 0.00000000 0.0000000 0.00000000 0.08302805 0.000000e+00
##  [6,] 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 1.308527e-01
##  [7,] 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 1.030822e-04
##  [8,] 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.000000e+00
##  [9,] 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 3.696276e-05
## [10,] 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.000000e+00
##               [,7]       [,8]         [,9]      [,10]
##  [1,] 0.0000000000 0.00000000 0.000000e+00 0.00000000
##  [2,] 0.0000000000 0.00000000 0.000000e+00 0.00000000
##  [3,] 0.0000000000 0.00000000 0.000000e+00 0.00000000
##  [4,] 0.0000000000 0.00000000 0.000000e+00 0.00000000
##  [5,] 0.0000000000 0.00000000 0.000000e+00 0.00000000
##  [6,] 0.0001030822 0.00000000 3.696276e-05 0.00000000
##  [7,] 0.1217530179 0.00000000 4.861388e-04 0.00000000
##  [8,] 0.0000000000 0.09427377 0.000000e+00 0.00000000
##  [9,] 0.0004861388 0.00000000 1.271970e-01 0.00000000
## [10,] 0.0000000000 0.00000000 0.000000e+00 0.09151835

a$wi[1:10,1:10]

##           [,1]     [,2]     [,3]    [,4]     [,5]     [,6]     [,7]    [,8]
##  [1,] 12.15457  0.00000 0.000000  0.0000  0.00000 0.000000 0.000000  0.0000
##  [2,]  0.00000 10.91906 0.000000  0.0000  0.00000 0.000000 0.000000  0.0000
##  [3,]  0.00000  0.00000 9.428381  0.0000  0.00000 0.000000 0.000000  0.0000
##  [4,]  0.00000  0.00000 0.000000 11.6498  0.00000 0.000000 0.000000  0.0000
##  [5,]  0.00000  0.00000 0.000000  0.0000 12.04412 0.000000 0.000000  0.0000
##  [6,]  0.00000  0.00000 0.000000  0.0000  0.00000 7.642636 0.000000  0.0000
##  [7,]  0.00000  0.00000 0.000000  0.0000  0.00000 0.000000 8.305836  0.0000
##  [8,]  0.00000  0.00000 0.000000  0.0000  0.00000 0.000000 0.000000 10.6074
##  [9,]  0.00000  0.00000 0.000000  0.0000  0.00000 0.000000 0.000000  0.0000
## [10,]  0.00000  0.00000 0.000000  0.0000  0.00000 0.000000 0.000000  0.0000
##           [,9]    [,10]
##  [1,] 0.000000  0.00000
##  [2,] 0.000000  0.00000
##  [3,] 0.000000  0.00000
##  [4,] 0.000000  0.00000
##  [5,] 0.000000  0.00000
##  [6,] 0.000000  0.00000
##  [7,] 0.000000  0.00000
##  [8,] 0.000000  0.00000
##  [9,] 7.873255  0.00000
## [10,] 0.000000 10.92677

## No. of detected edges
fit.gl=abs(a$wi)>1e-6
sum(fit.gl==1)/2

## [1] 205

## degree
glasso.degree=apply(fit.gl, 2, sum)

tab<- cbind.data.frame("Nodes"=colnames(data), "degree"=glasso.degree)
ss<- sort(tab$degree, ind=T, decreasing = T)
tab.sort<- tab[ss$ix,]
tab.sort$Nodes<- factor(tab.sort$Nodes, levels=unique(tab.sort$Nodes))
dat<- tab.sort[1:30,]

#ggplot(dat, aes(y=degree, x=Nodes)) +geom_bar(position="stack", stat="identity") + theme_bw() + theme(axis.text.x = element_text(colour="black",size=15,angle=45,hjust=1,vjust=1,face="plain"),axis.title.x = element_text(colour="black",size=25,angle=0,hjust=.5,vjust=.5,face="plain"), axis.text.y = element_text(colour="black",size=15,angle=0,hjust=1,vjust=0,face="plain"),  axis.title.y = element_text(colour="black",size=25,angle=90,hjust=.5,vjust=.5,face="plain"),legend.title = element_text(colour="black",size=15), legend.text=element_text(colour="black",size=15), legend.direction = "horizontal", legend.position = "top") + theme(plot.margin = unit(c(1,1,1,1), "cm"))

Running space.neighbor

load("data_retro.RData")  
data<- t(data)
p<- dim(data)[2] # no. of nodes
n=nrow(data)

alpha=1
l1=1/sqrt(n)*qnorm(1-alpha/(2*p^2))
iter=10

result1=space.neighbor(data, lam1=l1*0.7, lam2=0)
result1$ParCor[1:10,1:10]

##       [,1] [,2] [,3] [,4] [,5]       [,6] [,7] [,8]       [,9] [,10]
##  [1,]    1    0    0    0    0 0.00000000    0    0 0.00000000     0
##  [2,]    0    1    0    0    0 0.00000000    0    0 0.00000000     0
##  [3,]    0    0    1    0    0 0.00000000    0    0 0.00000000     0
##  [4,]    0    0    0    1    0 0.00000000    0    0 0.00000000     0
##  [5,]    0    0    0    0    1 0.00000000    0    0 0.00000000     0
##  [6,]    0    0    0    0    0 1.00000000    0    0 0.06109301     0
##  [7,]    0    0    0    0    0 0.00000000    1    0 0.00000000     0
##  [8,]    0    0    0    0    0 0.00000000    0    1 0.00000000     0
##  [9,]    0    0    0    0    0 0.06109301    0    0 1.00000000     0
## [10,]    0    0    0    0    0 0.00000000    0    0 0.00000000     1

## No. of detected edges
fit.space.nb=abs(result1$ParCor)>1e-6
sum(fit.space.nb==1)/2

## [1] 396

## degree
spacenb.degree=apply(fit.space.nb, 2, sum)

tab<- cbind.data.frame("Nodes"=colnames(data), "degree"=spacenb.degree)
ss<- sort(tab$degree, ind=T, decreasing = T)
tab.sort<- tab[ss$ix,]
tab.sort$Nodes<- factor(tab.sort$Nodes, levels=unique(tab.sort$Nodes))
dat<- tab.sort[1:30,]

#ggplot(dat, aes(y=degree, x=Nodes)) +geom_bar(position="stack", stat="identity") + theme_bw() + theme(axis.text.x = element_text(colour="black",size=15,angle=45,hjust=1,vjust=1,face="plain"),axis.title.x = element_text(colour="black",size=25,angle=0,hjust=.5,vjust=.5,face="plain"), axis.text.y = element_text(colour="black",size=15,angle=0,hjust=1,vjust=0,face="plain"),  axis.title.y = element_text(colour="black",size=25,angle=90,hjust=.5,vjust=.5,face="plain"),legend.title = element_text(colour="black",size=15), legend.text=element_text(colour="black",size=15), legend.direction = "horizontal", legend.position = "top") + theme(plot.margin = unit(c(1,1,1,1), "cm"))

Running space.neighbor

load("data_retro.RData")  
data<- t(data)
p<- dim(data)[2] # no. of nodes
n=nrow(data)

alpha=1
l1=1/sqrt(n)*qnorm(1-alpha/(2*p^2))
iter=10

result2=space.joint(data, lam1=6, lam2=0, iter=iter)

## [1] "iter=1"
## [1] "iter=2"
## [1] "iter=3"
## [1] "iter=4"
## [1] "iter=5"
## [1] "iter=6"
## [1] "iter=7"
## [1] "iter=8"
## [1] "iter=9"
## [1] "iter=10"

result2$ParCor[1:10,1:10]

##       [,1] [,2] [,3] [,4] [,5]       [,6] [,7] [,8]       [,9] [,10]
##  [1,]    1    0    0    0    0 0.00000000    0    0 0.00000000     0
##  [2,]    0    1    0    0    0 0.00000000    0    0 0.00000000     0
##  [3,]    0    0    1    0    0 0.00000000    0    0 0.00000000     0
##  [4,]    0    0    0    1    0 0.00000000    0    0 0.00000000     0
##  [5,]    0    0    0    0    1 0.00000000    0    0 0.00000000     0
##  [6,]    0    0    0    0    0 1.00000000    0    0 0.01794755     0
##  [7,]    0    0    0    0    0 0.00000000    1    0 0.00000000     0
##  [8,]    0    0    0    0    0 0.00000000    0    1 0.00000000     0
##  [9,]    0    0    0    0    0 0.01794755    0    0 1.00000000     0
## [10,]    0    0    0    0    0 0.00000000    0    0 0.00000000     1

## No. of detected edges
fit.space.jt=abs(result2$ParCor)>1e-6
sum(fit.space.jt==1)/2

## [1] 228

## degree
spacejt.degree=apply(fit.space.jt, 2, sum)

tab<- cbind.data.frame("Nodes"=colnames(data), "degree"=spacejt.degree)
ss<- sort(tab$degree, ind=T, decreasing = T)
tab.sort<- tab[ss$ix,]
tab.sort$Nodes<- factor(tab.sort$Nodes, levels=unique(tab.sort$Nodes))
dat<- tab.sort[1:30,]

#ggplot(dat, aes(y=degree, x=Nodes)) +geom_bar(position="stack", stat="identity") + theme_bw() + theme(axis.text.x = element_text(colour="black",size=15,angle=45,hjust=1,vjust=1,face="plain"),axis.title.x = element_text(colour="black",size=25,angle=0,hjust=.5,vjust=.5,face="plain"), axis.text.y = element_text(colour="black",size=15,angle=0,hjust=1,vjust=0,face="plain"),  axis.title.y = element_text(colour="black",size=25,angle=90,hjust=.5,vjust=.5,face="plain"),legend.title = element_text(colour="black",size=15), legend.text=element_text(colour="black",size=15), legend.direction = "horizontal", legend.position = "top") + theme(plot.margin = unit(c(1,1,1,1), "cm"))